home *** CD-ROM | disk | FTP | other *** search
/ SPACE 1 / SPACE - Library 1 - Volume 1.iso / program / 441 / fplib20 / log.c < prev    next >
C/C++ Source or Header  |  1990-11-23  |  1KB  |  46 lines
  1. /* LOG functions for Sozobon C.                        */
  2. /* Copyright © David Brooks, 1989 All Rights Reserved            */
  3.  
  4. #include <fplib.h>
  5. #include <errno.h>
  6.  
  7. /* LOG:                                    */
  8. /*                                    */
  9. /* This yanks out the exponent, performs a polynomial valid over +/-    */
  10. /* sqrt(2), and adds the scaled exponent back in.  Domain error on    */
  11. /* arguments <= 0.                            */
  12.  
  13. float log(a) fstruct a;
  14. {    register float    t1, t2, t3;
  15.     static fstruct nhuge = {0xFFFFFFFFL};    /* Largest negative    */
  16.  
  17.     if (a.f <= 0.0)
  18.     {    errno = EDOM;            /* Impossible really    */
  19.         return nhuge.f;
  20.     }
  21.  
  22.     if (a.f < M_SQRT2 && a.f > M_SQRT1_2)
  23.     {    t3 = 0.0;    /* arg between sqrt(2) and sqrt(2)/2    */
  24.         t2 = (t1 = a.f - 1.0) + 2.0;
  25.     }
  26.     else
  27.     {    t3 = (float)(((int)a.sc[3] << 1) - 0x81) * 0.34657359;
  28.             /* ...(log2(a) * 2 + 1) * ln(2)/2 */
  29.         a.sc[3] = BIAS;            /* Scale back to range */
  30.         t2 = a.f + M_SQRT1_2;
  31.         t1 = a.f - M_SQRT1_2;
  32.     }
  33.  
  34.     t1 /= t2;
  35.     t2 = t1 * t1;
  36.  
  37.     return t3 + t1 * (2.0000008 + t2 * (0.66644078 + t2 * 0.41517739));
  38. }
  39.  
  40. /* LOG10 is easy. */
  41.  
  42. float log10(a) float a;
  43. {
  44.     return M_LOG10E * log(a);
  45. }
  46.